Hands-on Exercise 3 and 4

Published

January 19, 2024

Introduction to Part 1: 1st Order Spatial Point Patterns Analysis Methods

In this exercse, we will be learning about 1st order Spatial Point Patterns Analysis

Setting Up

The Data

The following datasets will be used

  • CHILDCARE which was downloaded from data.gov.sg
  • 2014 Master Plan Subzone Boundary (No Sea) which was downloaded from data.gov.sg
  • CoastalOutline, which was derived from 2019 Master Plan Subzone Boundary (No Sea). It was downloaded from data.gov.sg

Downloading Packages

The following packages will be used:

  • maptools
  • sf
  • raster
  • spatstat
  • tmap
Show the code
pacman::p_load(maptools, sf, raster, spatstat, tmap)

As there may be issues with installing maptools, you may use the following:

Show the code
install.packages("maptools", repos = "https://packagemanager.posit.co/cran/2023-10-13")

Spatial Data Wrangling

Importing the Spatial Data

Import the spatial data using st_read().

Show the code
childcare_sf <- st_read("data/child-care-services-geojson.geojson") %>%
  st_transform(crs = 3414)
Reading layer `child-care-services-geojson' from data source 
  `D:\KrisLBT\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data\child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
Show the code
mpsz_sf <- st_read(dsn = "data", 
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `D:\KrisLBT\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
Show the code
sg_sf <-
  st_read(dsn="data",
          layer= "CostalOutline")
Reading layer `CostalOutline' from data source 
  `D:\KrisLBT\IS415-GAA\Hands-on_Ex\Hands-on_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21

DIY

From here, we will check the referencing system.

Show the code
st_crs(childcare_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
Show the code
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Show the code
st_crs(sg_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

childcare_sf and mpsz_sf both do not have the correct crs so we will need to set it to CRS = 3414.

DIY 2

Now, I will be setting the CRS to the correct one.Both are in SVY21/Singapore TM

Show the code
mpsz_sf <- st_set_crs(mpsz_sf, 3414)
Show the code
sg_sf <- st_set_crs(sg_sf, 3414)

Mapping the geospatial data sets

We can plot a map to show the spatial patterns of our data.

Show the code
tm_shape(mpsz_sf)+
  tm_polygons() +
tm_shape(childcare_sf) +
  tm_dots()

We can also create a pin map as shown below:

Show the code
tmap_mode("view")
tm_shape(childcare_sf)+
  tm_dots()
Show the code
tmap_mode("plot")
Caution

Always remember to switch back to plot mode after the interactive map. This is because, each interactive mode will consume a connection. You should also avoid displaying ecessive numbers of interactive maps (i.e. not more than 10) in one RMarkdown document when publishing on Netlify.

Geospatial Data wrangling

Although simple feature data frame is gaining popularity again sp’s Spatial* classes, there are, however, many geospatial analysis packages require the input geospatial data in sp’s Spatial* classes. In this section, you will learn how to convert simple feature data frame to sp’s Spatial* class.

Converting sf data frames to sp’s Spatial* class –> idk

Show the code
childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)
Show the code
print(mpsz)
class       : SpatialPolygonsDataFrame 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 15
names       : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C,          INC_CRC, FMEL_UPD_D,     X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
min values  :        1,          1, ADMIRALTY,    AMSZ01,      N, ANG MO KIO,         AM, CENTRAL REGION,       CR, 00F5E30B5C9B7AD8,      16409,  5092.8949,  19579.069, 871.554887798, 39437.9352703 
max values  :      323,         17,    YUNNAN,    YSSZ09,      Y,     YISHUN,         YS,    WEST REGION,       WR, FFCCF172717C2EAF,      16409, 50424.7923, 49552.7904, 68083.9364708,  69748298.792 
Show the code
print(childcare)
class       : SpatialPointsDataFrame 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>018989</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>1, MARINA BOULEVARD, #B1 - 01, ONE MARINA BOULEVARD, SINGAPORE 018989</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>THE LITTLE SKOOL-HOUSE INTERNATIONAL PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>08F73931F4A691F4</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                  <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>829646</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>200, PONGGOL SEVENTEENTH AVENUE, SINGAPORE 829646</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>RAFFLES KIDZ @ PUNGGOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>379D017BF244B0FA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
Show the code
print(sg)
class       : SpatialPolygonsDataFrame 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       : GDO_GID, MSLINK, MAPID,              COSTAL_NAM 
min values  :       1,      1,     0,             ISLAND LINK 
max values  :      60,     67,     0, SINGAPORE - MAIN ISLAND 

Converting the Spatial* class into generic sp format

Show the code
childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")
Show the code
print(childcare_sp)
class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
Show the code
print(sg_sp)
class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Converting the Spatial* class into generic sp format

The code chunk below converts the Spatial* classes into generic sp objects:

Show the code
childcare_sp <- as(childcare, "SpatialPoints")
sg_sg <- as(sg, "SpatialPolygons")

Next, I will display the sp objects as shown below:

Show the code
childcare_sp
class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
Show the code
sg_sp
class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Answer to challenge: > The class Spatial only holds metadata common to all derived classes (bounding box,coordinate reference system),and is convenient for defining methods that are common to all derived classes.

– From here

Converting the generic sp format into spatstat’s ppp format

About ppp format

Conversion:

Show the code
childcare_ppp<- as(childcare_sp,"ppp")
childcare_ppp
Planar point pattern: 1545 points
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units

Now, I will plot childcare_ppp and examine the difference:

Show the code
plot(childcare_ppp)

Summary statistics:

Show the code
summary(childcare_ppp)
Planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units

Handling duplicated points

We can check for duplication in a ppp object by using the code chunk below:

Show the code
any(duplicated(childcare_ppp))
[1] TRUE

To count the number of co-incidence point, we will use the multiplicity() function as shown below (as there are too many entries, I have set eval to FALSE):

Show the code
multiplicity(childcare_ppp)

To find out how many locations have more than 1 point,event, we can use the code chunk below:

Show the code
sum(multiplicity(childcare_ppp)>1)
[1] 128

This shwos that there are 1128duplicated point events

To view the locations of these duplicate point events, we will plot childcare data by using the following code chunk:

Show the code
tmap_mode("view")
tm_shape(childcare) +
  tm_dots(alpha=0.4,
          size = 0.05)
Show the code
tmap_mode('plot')

Answer to challenge > The ones with multiplicity are the ones that are darker compared to others

There are three ways to solve the problem: - Delete the duplicates (however, you may lose some useful events) - jittering, which will add a small perturbation to the duplicate points so they do not occupy the exact same space. - Make each point “unique” and the attach the duplicates of the points to the patterns as marks, as attributes of the points.

The code chunk below shows the jittering approach (method 3)

Show the code
childcare_ppp_jit <- rjitter(childcare_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)

Now, I will check if there remains any duplicated data

Show the code
any(duplicated(childcare_ppp_jit))
[1] FALSE

Alternatively,

Show the code
sum(multiplicity(childcare_ppp_jit)>1)
[1] 0

Creating own owin object

When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.

Show the code
sg_owin <- as(sg_sp, "owin")

the output object can be displayed using the plot() functions:

Show the code
plot(sg_owin)

And here is the summary:

Show the code
summary(sg_owin)
Window: polygonal boundary
60 separate polygons (no holes)
            vertices        area relative.area
polygon 1         38 1.56140e+04      2.09e-05
polygon 2        735 4.69093e+06      6.27e-03
polygon 3         49 1.66986e+04      2.23e-05
polygon 4         76 3.12332e+05      4.17e-04
polygon 5       5141 6.36179e+08      8.50e-01
polygon 6         42 5.58317e+04      7.46e-05
polygon 7         67 1.31354e+06      1.75e-03
polygon 8         15 4.46420e+03      5.96e-06
polygon 9         14 5.46674e+03      7.30e-06
polygon 10        37 5.26194e+03      7.03e-06
polygon 11        53 3.44003e+04      4.59e-05
polygon 12        74 5.82234e+04      7.78e-05
polygon 13        69 5.63134e+04      7.52e-05
polygon 14       143 1.45139e+05      1.94e-04
polygon 15       165 3.38736e+05      4.52e-04
polygon 16       130 9.40465e+04      1.26e-04
polygon 17        19 1.80977e+03      2.42e-06
polygon 18        16 2.01046e+03      2.69e-06
polygon 19        93 4.30642e+05      5.75e-04
polygon 20        90 4.15092e+05      5.54e-04
polygon 21       721 1.92795e+06      2.57e-03
polygon 22       330 1.11896e+06      1.49e-03
polygon 23       115 9.28394e+05      1.24e-03
polygon 24        37 1.01705e+04      1.36e-05
polygon 25        25 1.66227e+04      2.22e-05
polygon 26        10 2.14507e+03      2.86e-06
polygon 27       190 2.02489e+05      2.70e-04
polygon 28       175 9.25904e+05      1.24e-03
polygon 29      1993 9.99217e+06      1.33e-02
polygon 30        38 2.42492e+04      3.24e-05
polygon 31        24 6.35239e+03      8.48e-06
polygon 32        53 6.35791e+05      8.49e-04
polygon 33        41 1.60161e+04      2.14e-05
polygon 34        22 2.54368e+03      3.40e-06
polygon 35        30 1.08382e+04      1.45e-05
polygon 36       327 2.16921e+06      2.90e-03
polygon 37       111 6.62927e+05      8.85e-04
polygon 38        90 1.15991e+05      1.55e-04
polygon 39        98 6.26829e+04      8.37e-05
polygon 40       415 3.25384e+06      4.35e-03
polygon 41       222 1.51142e+06      2.02e-03
polygon 42       107 6.33039e+05      8.45e-04
polygon 43         7 2.48299e+03      3.32e-06
polygon 44        17 3.28303e+04      4.38e-05
polygon 45        26 8.34758e+03      1.11e-05
polygon 46       177 4.67446e+05      6.24e-04
polygon 47        16 3.19460e+03      4.27e-06
polygon 48        15 4.87296e+03      6.51e-06
polygon 49        66 1.61841e+04      2.16e-05
polygon 50       149 5.63430e+06      7.53e-03
polygon 51       609 2.62570e+07      3.51e-02
polygon 52         8 7.82256e+03      1.04e-05
polygon 53       976 2.33447e+07      3.12e-02
polygon 54        55 8.25379e+04      1.10e-04
polygon 55       976 2.33447e+07      3.12e-02
polygon 56        61 3.33449e+05      4.45e-04
polygon 57         6 1.68410e+04      2.25e-05
polygon 58         4 9.45963e+03      1.26e-05
polygon 59        46 6.99702e+05      9.35e-04
polygon 60        13 7.00873e+04      9.36e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 748741000 square units
Fraction of frame area: 0.414

Combining point events object and owin object

In the last step of geospatial wrangling, we will extract childcare events that are located within Singapore

Show the code
childcareSG_ppp = childcare_ppp[sg_owin]

DIY

Show the code
plot(childcareSG_ppp)

First-order Spatial Point Patterns Analysis

In this section, I will perform first-order SPPA by using spatstat package. The hands-on exercise will focus on:

  • deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes,
  • performing Confirmatory Spatial Point Patterns Analysis by using Nearest Neighbour statistics.

Kernel Density Estimation

Computing kernel density estimation using automatic bandwidth selection method

The code chunk below computes a kernel density by using the following configurations of density() of spatstat:

  • bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl().
  • The smoothing kernel used is gaussian, which is the default. Other smoothing methods are: “epanechnikov”, “quartic” or “disc”.
  • The intensity estimate is corrected for edge effect bias by using method described by Jones (1993) and Diggle (2010, equation 18.9). The default is FALSE.
Show the code
kde_childcareSG_bw <- density(childcareSG_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

The below code chunk will display the derived kernel display.

Show the code
plot(kde_childcareSG_bw)

The density values of the output range from 0 to 0.000035 which is way too small to comprehend. This is because the default unit of measurement of svy21 is in meter. As a result, the density values computed is in “number of points per square meter”.

Before we move on to next section, it is good to know that you can retrieve the bandwidth used to compute the kde layer by using the code chunk below.

Show the code
bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
298.4095 

Rescaling KDE values

In the code chunk below, rescale() is used to covert the unit of measurement from meter to kilometer.

Show the code
childcareSG_ppp.km <- rescale(childcareSG_ppp, 1000, "km")

We can re-run density() using the rescaled dataset and plot the output kde map.

Show the code
kde_childcareSG.bw <- density(childcareSG_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG.bw)

Notice that output image looks identical to the earlier version, the only changes in the data values (refer to the legend).

Working with different automatic badwidth methods

Beside bw.diggle(), there are three other spatstat functions can be used to determine the bandwidth, they are: bw.CvL(), bw.scott(), and bw.ppl().

Let us take a look at the bandwidth return by these automatic bandwidth calculation methods by using the code chunk below.

Show the code
bw.CvL(childcareSG_ppp.km)
   sigma 
4.543278 
Show the code
bw.scott(childcareSG_ppp.km)
 sigma.x  sigma.y 
2.224898 1.450966 
Show the code
bw.ppl(childcareSG_ppp.km)
    sigma 
0.3897114 
Show the code
bw.diggle(childcareSG_ppp.km)
    sigma 
0.2984095 

Baddeley et. (2016) suggested the use of the bw.ppl() algorithm because in ther experience it tends to produce the more appropriate values when the pattern consists predominantly of tight clusters. But they also insist that if the purpose of once study is to detect a single tight cluster in the midst of random noise then the bw.diggle() method seems to work best.

The code chunk beow will be used to compare the output of using bw.diggle and bw.ppl methods.

Show the code
kde_childcareSG.ppl <- density(childcareSG_ppp.km, 
                               sigma=bw.ppl, 
                               edge=TRUE,
                               kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

Working with different kernel methods

By default, the kernel method used in density.ppp() is gaussian. But there are three other options, namely: Epanechnikov, Quartic and Dics.

The code chunk below will be used to compute three more kernel density estimations by using these three kernel function.

Show the code
par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="gaussian"), 
     main="Gaussian")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="epanechnikov"), 
     main="Epanechnikov")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="quartic"), 
     main="Quartic")
plot(density(childcareSG_ppp.km, 
             sigma=bw.ppl, 
             edge=TRUE, 
             kernel="disc"), 
     main="Disc")

Fixed and Adaptive KDE

Computing KDE by using fixed bandwidth

Next, you will compute a KDE layer by defining a bandwidth of 600 meter. Notice that in the code chunk below, the sigma value used is 0.6. This is because the unit of measurement of childcareSG_ppp.km object is in kilometer, hence the 600m is 0.6km.

Show the code
kde_childcareSG_600 <- density(childcareSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600)

Computing KDE by using adaptive bandwidth

Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units for example urban versus rural. One way to overcome this problem is by using adaptive bandwidth instead.

In this section, you will learn how to derive adaptive kernel density estimation by using density.adaptive() of spatstat.

Show the code
kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)

We can compare the fixed and adaptive kernel density estimation outputs by using the code chunk below.

Show the code
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")

Converting KDE output into grid object.

The result is the same, we just convert it so that it is suitable for mapping purposes

Show the code
gridded_kde_childcareSG_bw <- as.SpatialGridDataFrame.im(kde_childcareSG.bw)
spplot(gridded_kde_childcareSG_bw)

Converting gridded output into raster

Next, we will convert the gridded kernal density objects into RasterLayer object by using raster() of raster package.

Show the code
kde_childcareSG_bw_raster <- raster(gridded_kde_childcareSG_bw)

Let us take a look at the properties:

Show the code
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -8.476185e-15, 28.51831  (min, max)

Notice that the crs property is NA.

Assigning projections system

The code chunk below will be used to include the CRS information on kde_childcareSG_bw_raster Raster Layer.

Show the code
projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -8.476185e-15, 28.51831  (min, max)

Visualising the output in tmap

Finally, we will display the raster in cartographic quality map using tmap package.

Show the code
tm_shape(kde_childcareSG_bw_raster) +
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom", frame = FALSE))

Notice that the raster values are encoded explicitly onto the raster pixel using the values in “v”” field.

Comparing Spatial Point Patterns using KDE

In this section, you will learn how to compare KDE of childcare at Ponggol, Tampines, Chua Chu Kang and Jurong West planning areas.

Extracting study area

The code chunk below will be used to extract the target planning areas.

Show the code
pg = mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]
tm = mpsz[mpsz@data$PLN_AREA_N == "TAMPINES",]
ck = mpsz[mpsz@data$PLN_AREA_N == "CHOA CHU KANG",]
jw = mpsz[mpsz@data$PLN_AREA_N == "JURONG WEST",]

Plotting target planning areas

Show the code
par(mfrow=c(2,2))
plot(pg, main = "Ponggol")
plot(tm, main = "Tampines")
plot(ck, main = "Choa Chu Kang")
plot(jw, main = "Jurong West")

Converting the spatial point data frame into generic sp format

Next, we will convert these SpatialPolygonsDataFrame layers into generic spatialpolygons layers.

Show the code
pg_sp = as(pg, "SpatialPolygons")
tm_sp = as(tm, "SpatialPolygons")
ck_sp = as(ck, "SpatialPolygons")
jw_sp = as(jw, "SpatialPolygons")

Creating owin object

Now, we will convert these SpatialPolygons objects into owin objects that is required by spatstat.

Show the code
pg_owin = as(pg_sp, "owin")
tm_owin = as(tm_sp, "owin")
ck_owin = as(ck_sp, "owin")
jw_owin = as(jw_sp, "owin")

Combining childcare points and the study area

By using the code chunk below, we are able to extract childcare that is within the specific region to do our analysis later on.

Show the code
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]

Next, rescale() function is used to trasnform the unit of measurement from metre to kilometre.

Show the code
childcare_pg_ppp.km = rescale(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale(childcare_jw_ppp, 1000, "km")

The code chunk below is used to plot these four study areas and the locations of the childcare centres.

Show the code
par(mfrow=c(2,2))
plot(childcare_pg_ppp.km, main="Punggol")
plot(childcare_tm_ppp.km, main="Tampines")
plot(childcare_ck_ppp.km, main="Choa Chu Kang")
plot(childcare_jw_ppp.km, main="Jurong West")

Computing KDE

The code chunk below will be used to compute the KDE of these four planning area. bw.diggle method is used to derive the bandwidth of each

Show the code
par(mfrow=c(2,2))
plot(density(childcare_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tempines")
plot(density(childcare_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="JUrong West")

Computing fixed bandwidth KDE

For comparison purposes, we will use 250m as the bandwidth.

Show the code
par(mfrow=c(2,2))
plot(density(childcare_ck_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="JUrong West")
plot(density(childcare_pg_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km, 
             sigma=0.25, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")

Nearest Neighbour Analysis

In this section, we will perform the Clark-Evans test of aggregation for a spatial point pattern by using clarkevans.test() of statspat.

The test hypotheses are:

Ho = The distribution of childcare services are randomly distributed.

H1= The distribution of childcare services are not randomly distributed.

The 95% confident interval will be used.

Testing spatial point patterns using Clark and Evans Test

Show the code
clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.54756, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

Conclusion: “We reject the null hypothesis”

Clark and Evans Test: Choa Chu Kang planning area

In the code chunk below, clarkevans.test() of spatstat is used to performs Clark-Evans test of aggregation for childcare centre in Choa Chu Kang planning area.

Show the code
clarkevans.test(childcare_ck_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.95209, p-value = 0.4741
alternative hypothesis: two-sided

Conclusion: we do not reject null hypothesis

Clark and Evans Test: Tampines planning area

Show the code
clarkevans.test(childcare_tm_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.78152, p-value = 8.042e-05
alternative hypothesis: two-sided

Conclusion: We reject the null hypothesis

Introduction to Part 2: 2nd Order Spatial Point Patterns Analysis

Analysiing Spatial Point Process Using G-Function

The G function measures the distribution of the distances from an arbitrary event to its nearest event. In this section, you will learn how to compute G-function estimation by using Gest() of spatstat package. You will also learn how to perform monta carlo simulation test using envelope() of spatstat package.

Choa Chu Kang planning area

Computing G-function estimation

The code chunk below is used to compute G-function using Gest() of spatat package.

Show the code
G_CK = Gest(childcare_ck_ppp, correction = "border")
plot(G_CK, xlim=c(0,500))

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Monte Carlo test with G-fucntion

Show the code
G_CK.csr <- envelope(childcare_ck_ppp, Gest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
Show the code
plot(G_CK.csr)

Tampines planning area

Computing G-function estimation

Show the code
G_tm = Gest(childcare_tm_ppp, correction = "best")
plot(G_tm)

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Tampines are randomly distributed.

H1= The distribution of childcare services at Tampines are not randomly distributed.

The null hypothesis will be rejected is p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

Show the code
G_tm.csr <- envelope(childcare_tm_ppp, Gest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
Show the code
plot(G_tm.csr)

Analysing Spatial Point Process Using F-Function

The F function estimates the empty space function F(r) or its hazard rate h(r) from a point pattern in a window of arbitrary shape. In this section, you will learn how to compute F-function estimation by using Fest() of spatstat package. You will also learn how to perform monta carlo simulation test using envelope() of spatstat package.

Choa Chu Kang planning area

Computing F-function estimation

The code chunk below is used to compute F-function using Fest() of spatat package.

Show the code
F_CK = Fest(childcare_ck_ppp)
plot(F_CK)

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Monte Carlo test with F-fucntion

Show the code
F_CK.csr <- envelope(childcare_ck_ppp, Fest, nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
Show the code
plot(F_CK.csr)

Tampines planning area

Computing F-function estimation

Monte Carlo test with F-function

Show the code
F_tm = Fest(childcare_tm_ppp, correction = "best")
plot(F_tm)

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Tampines are randomly distributed.

H1= The distribution of childcare services at Tampines are not randomly distributed.

The null hypothesis will be rejected is p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

Show the code
F_tm.csr <- envelope(childcare_tm_ppp, Fest, correction = "all", nsim = 999)
Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60..
.......70.........80.........90.........100.........110.........120.........130
.........140.........150.........160.........170.........180.........190........
.200.........210.........220.........230.........240.........250.........260......
...270.........280.........290.........300.........310.........320.........330....
.....340.........350.........360.........370.........380.........390.........400..
.......410.........420.........430.........440.........450.........460.........470
.........480.........490.........500.........510.........520.........530........
.540.........550.........560.........570.........580.........590.........600......
...610.........620.........630.........640.........650.........660.........670....
.....680.........690.........700.........710.........720.........730.........740..
.......750.........760.........770.........780.........790.........800.........810
.........820.........830.........840.........850.........860.........870........
.880.........890.........900.........910.........920.........930.........940......
...950.........960.........970.........980.........990........
999.

Done.
Show the code
plot(F_tm.csr)

Analysing Spatial Point Process Using K-Function

K-function measures the number of events found up to a given distance of any particular event. In this section, you will learn how to compute K-function estimates by using Kest() of spatstat package. You will also learn how to perform monta carlo simulation test using envelope() of spatstat package.

Choa Chu Kang planning area

Computing K-function estimate

Show the code
K_ck = Kest(childcare_ck_ppp, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(m)")

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

Show the code
K_ck.csr <- envelope(childcare_ck_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
Show the code
plot(K_ck.csr, . - r ~ r, xlab="d", ylab="K(d)-r")

Tampines planning area

Computing K-function estimation

Show the code
K_tm = Kest(childcare_tm_ppp, correction = "Ripley")
plot(K_tm, . -r ~ r, 
     ylab= "K(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Tampines are randomly distributed.

H1= The distribution of childcare services at Tampines are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

Show the code
K_tm.csr <- envelope(childcare_tm_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
Show the code
plot(K_tm.csr, . - r ~ r, 
     xlab="d", ylab="K(d)-r", xlim=c(0,500))

Analysing Spatial Point Process Using L-Function

In this section, you will learn how to compute L-function estimation by using Lest() of spatstat package. You will also learn how to perform monta carlo simulation test using envelope() of spatstat package.

Choa Chu Kang planning area

Show the code
L_ck = Lest(childcare_ck_ppp, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Choa Chu Kang are randomly distributed.

H1= The distribution of childcare services at Choa Chu Kang are not randomly distributed.

The null hypothesis will be rejected if p-value if smaller than alpha value of 0.001.

The code chunk below is used to perform the hypothesis testing.

Show the code
L_ck.csr <- envelope(childcare_ck_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
Show the code
plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

Tampines planning area

Computing L-function estimate

Show the code
L_tm = Lest(childcare_tm_ppp, correction = "Ripley")
plot(L_tm, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)", 
     xlim=c(0,1000))

Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at Tampines are randomly distributed.

H1= The distribution of childcare services at Tampines are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

The code chunk below will be used to perform the hypothesis testing.

Show the code
L_tm.csr <- envelope(childcare_tm_ppp, Lest, nsim = 99, rank = 1, glocal=TRUE)
Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.
Show the code
plot(L_tm.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500))